# TODO: Add comment
# 
# Author: guochun
###############################################################################


oneSampling.tSquare=function(distanceOBJ,populationData){
	#get a nearest event from a random point
	availible=randomPointToEvent(distanceOBJ,populationData)
	d=as.numeric(attr(availible,"d"))
	#calculate the qualified area
	if(d[4]!=d[2]){
		tgbeta=(d[1]-d[3])/(d[4]-d[2])
		b=d[4]-tgbeta*d[3]
		yo=selectionLine(d[1],tgbeta,b)
		sLy=selectionLine(populationData@x,tgbeta,b)
		if(d[2]<yo){
			selection=which(populationData@y>sLy)
		}else{
			selection=which(populationData@y<sLy)
		}
	}else{
		#if the first selected two point is on the sample horizental line
		if(d[1]>d[3]){
			selection=which(populationData@x<d[3])
		}else{
			selection=which(populationData@x>d[3])
		} 
	}
	#get the second nearest event form the previous event from the qualified area
	window=owin(c(0,populationData@plotdim[1]),c(0,populationData@plotdim[2]))   
	e1=ppp(x=d[3],y=d[4],window=window,check=F)
	e2=ppp(x=populationData@x[selection],y=populationData@y[selection],window=window,check=F)
	del=which(e2$x==d[3] & e2$y==d[4])
	if(length(del)!=0){
		e2=e2[-del]
	}
	nd=nncross(e1,e2)
	if(is.na(nd[1,2])){
		attr(availible,"d")=NA
	}else{
		attr(availible,"d")=data.frame(xfrom=d[1],yfrom=d[2],xBy=d[3],yBy=d[4],xTo=e2[nd[1,2]]$x,yTo=e2[nd[1,2]]$y,d1=d[5],d2=nd[1,1])	
	}
	
	return(availible)
}


setSampleResult.tSquare=function(sampleOBJ,selection,populationData){
	return(appendResult(sampleOBJ,selection,populationData))
}

samplingQuick.tSquare=function(sampleMethod,populationData){
	#browser()
	#多取20%的点，最后在删除重复点后，若有多余，则可以随机选择需要的点数
	if(sampleMethod@composite){
		if(sampleMethod@time%%sampleMethod@k != 0)
			stop("the total sample time can not be divided by 2")
		spt=round(sampleMethod@time/sampleMethod@k*1.2)
	}else{
		spt=round(sampleMethod@time*1.2)
	}
	if(spt>=length(populationData@x))
		spt=length(populationData@x)
	
	repl=TRUE
	rep.time=1
	while(repl & rep.time < sampleMethod@stopLimit){
		#define the row data matrix
		D=matrix(c(populationData@x,populationData@y),nrow=length(populationData@x))
		#random choose all points
		x=runif(spt, sampleMethod@edge, populationData@plotdim[1]-sampleMethod@edge)
		y=runif(spt, sampleMethod@edge, populationData@plotdim[2]-sampleMethod@edge)
		#define the first data matrix
		A=matrix(c(x,y),nrow=length(x))
		
		ki=knnx.index(D,A,k=1)[,1]
		kineg=which(ki == -1)
		if(length(kineg)>0){
			ki=ki[-kineg]
			A=A[-kineg,]
		}
			
		#define the second data matrix
		B=D[ki,]
		#modifiy the original D matrix
	    #不需要修改原始点了，因为这样会有误差，人为的稀疏了点，但为了避便找到与自己相同的点，需要调整k
		#D=D[-ki,]
		k=2
		final.C=rep(NA,dim(B)[1])
		needadd=rep(0,dim(B)[1])
		redo=TRUE
		
		kmax=(sampleMethod@k+1)*10
		#如果取临近个体到第20个还未找到符合要求的，则退出while
		kjtotal=knnx.index(D,B,k=kmax)
		while(any(redo) & k<kmax){
			kj=kjtotal[,k]
			#有些树可能会找不到临体，因此需要再考虑这个
			den=which(kj==-1)
			if(length(den)>0){
				kj[den]=1
			}
			C=D[kj,]
			
			angle=angleBAC(A[,1],A[,2],B[,1],B[,2],C[,1],C[,2])
			redo = angle <= pi/2
			if(any(is.na(redo)))
				browser()
	        redo[den]=TRUE
			needadd[!redo]=needadd[!redo]+1
			#这里表示只有到第k个邻体时，才要值
			add= needadd >=sampleMethod@k & is.na(final.C)
			
			final.C[add]=kj[add]
			redo = is.na(final.C)
			k=k+1
		}
		#先删除重复和NA的点
	    nNA=which(is.na(final.C))
		temp=table(final.C)
		#这里match时，只出现第一次出现的重复点
		nRep=match(as.numeric(names(temp)[which(temp>1)]),final.C)
		del=unique(c(nNA,nRep))
		nleft=spt-length(del)
		if(nleft>= (sampleMethod@time/sampleMethod@k)){
			rand.sl=sample(1:nleft,sampleMethod@time/sampleMethod@k)
			repl=FALSE
			if(length(del)>0){
				final.C=final.C[-del]
				B=B[-del,]
				if(sampleMethod@composite){
					A=A[-del,]
				}
			}
			
			final.C=final.C[rand.sl]
			B=B[rand.sl,]
			if(sampleMethod@composite){
				A=A[rand.sl,]
			}
		}else{
			#如果剩下的还不足够符合条件，则在重做一次
			repl=TRUE
			rep.time=rep.time+1
		}
		
		
	}
	if(rep.time < sampleMethod@stopLimit){
		final.C=D[final.C,]
		d=sqrt((B[,1]-final.C[,1])^2+(B[,2]-final.C[,2])^2)
		if(sampleMethod@composite){
			d2=sqrt((B[,1]-A[,1])^2+(B[,2]-A[,2])^2)
			d=c(d,d2)
		}
		
		return(d)
	}else{
		return(NA)
	}
	
}




# the angle of bac
angleBAC=function(x1,y1,x2,y2,x3,y3){
	b2=(x1-x2)^2+(y1-y2)^2
	c2=(x3-x2)^2+(y3-y2)^2
	a2=(x1-x3)^2+(y1-y3)^2
	bac=acos(round((b2+c2-a2)/(2*sqrt(b2*c2)),10))
	return(bac)
}
